Simple ellipsoid

import dolfin
from fenics_plotly import plot

import pulse

try:
    from dolfin_adjoint import Constant, DirichletBC, Function, Mesh, interpolate
except ImportError:
    from dolfin import Function, Constant, DirichletBC, Mesh, interpolate

gamma_space = "R_0"

geometry = pulse.HeartGeometry.from_file(pulse.mesh_paths["simple_ellipsoid"])
# geometry = pulse.geometries.prolate_ellipsoid_geometry(mesh_size_factor=3.0)

2023-10-15 20:16:28,140 [1193] INFO     pulse.geometry_utils: 
Load mesh from h5

if gamma_space == "regional":
    activation = pulse.RegionalParameter(geometry.cfun)
    target_activation = pulse.dolfin_utils.get_constant(0.2, len(activation))
else:
    activation = Function(dolfin.FunctionSpace(geometry.mesh, "R", 0))
    target_activation = Constant(0.2)

matparams = pulse.HolzapfelOgden.default_parameters()
material = pulse.HolzapfelOgden(
    activation=activation,
    parameters=matparams,
    f0=geometry.f0,
    s0=geometry.s0,
    n0=geometry.n0,
)

# LV Pressure
lvp = Constant(0.0)
lv_marker = geometry.markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]

# Add spring term at the epicardium of stiffness 1.0 kPa/cm^2 to represent pericardium
base_spring = 1.0
robin_bc = [
    pulse.RobinBC(value=Constant(base_spring), marker=geometry.markers["EPI"][0]),
]

# Fix the basal plane in the longitudinal direction
# 0 in V.sub(0) refers to x-direction, which is the longitudinal direction
def fix_basal_plane(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    bc = DirichletBC(
        V.sub(0),
        Constant(0.0),
        geometry.ffun,
        geometry.markers["BASE"][0],
    )
    return bc

dirichlet_bc = (fix_basal_plane,)

# Collect boundary conditions
bcs = pulse.BoundaryConditions(
    dirichlet=dirichlet_bc,
    neumann=neumann_bc,
    robin=robin_bc,
)

# Create the problem
problem = pulse.MechanicsProblem(geometry, material, bcs)

# Solve the problem
pulse.iterate.iterate(problem, (lvp, activation), (1.0, target_activation))

2023-10-15 20:16:28,353 [1193] INFO     pulse.iterate: Iterating to target control (f_36)...

2023-10-15 20:16:28,354 [1193] INFO     pulse.iterate: Current control: f_36 = 0.000,0.000

2023-10-15 20:16:28,355 [1193] INFO     pulse.iterate: Target: 1.000,0.200

2023-10-15 20:16:28,509 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,510 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,511 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,511 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,512 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,513 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,514 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,514 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,515 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,515 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:28,516 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:16:29,779 [1193] DEBUG    UFL_LEGACY: Blocks of each mode: 
  99	full

2023-10-15 20:16:30,049 [1193] DEBUG    UFL_LEGACY: Blocks of each mode: 
  18	full

2023-10-15 20:16:30,145 [1193] DEBUG    UFL_LEGACY: Blocks of each mode: 

2023-10-15 20:17:59,701 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,702 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,703 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,704 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,704 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,705 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,706 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,706 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,707 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,707 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,708 [1193] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.

2023-10-15 20:17:59,924 [1193] DEBUG    UFL_LEGACY: Blocks of each mode: 
  10	full

2023-10-15 20:18:00,176 [1193] DEBUG    UFL_LEGACY: Blocks of each mode: 
  3	full

2023-10-15 20:18:00,271 [1193] DEBUG    UFL_LEGACY: Blocks of each mode: 
  3	full

*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.

*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.

*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.

*** Warning: PETSc SNES solver diverged in 3 iterations with divergence reason DIVERGED_DTOL.

*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.

*** Warning: PETSc SNES solver diverged in 3 iterations with divergence reason DIVERGED_DTOL.

*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.

*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.

([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 57),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 300),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 375),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 442),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 509),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 568),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 627),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 726),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 785),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 868),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 991),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1098),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1205),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1280),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1332)],
 [(Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 53),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 55)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 296),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 298)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 371),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 373)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 438),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 440)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 505),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 507)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 564),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 566)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 623),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 625)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 722),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 724)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 781),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 783)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 864),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 866)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 987),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 989)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1094),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1096)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1201),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1203)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1276),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1278)),
  (Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 1328),
   Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1330))])

# Get the solution
u, p = problem.state.split(deepcopy=True)

# Move mesh accoring to displacement
u_int = interpolate(u, dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1))
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)

fig = plot(geometry.mesh, opacity=0.0, show=False)
fig.add_plot(plot(mesh, color="red", show=False))
fig.show()